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Abstract 


The solution of partial differential equations (PDEs) with Walsh functions offers new opportunities 
to simulate many challenging problems in mathematical physics. The approach was developed to 
better simulate hypersonic flows with shocks on unstructured grids. It is unique in that integrals and 
derivatives are computed using simple matrix multiplication of series representations of functions 
without the need for divided differences. The product of any two Walsh functions is another Walsh 
function - a feature that radically changes an algorithm for solving PDEs. A FORTRAN module 
for supporting Walsh function simulations is documented. A FORTRAN code is also documented 
with options for solving time-dependent problems: an advection equation, a Burgers equation, and 
a Riemann problem. The sample problems demonstrate the usage of the Walsh function module 
including such features as operator overloading, Fast Walsh Transforms in multi-dimensions, and a 
Fast Walsh reciprocal. 
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1 Introduction 


The solution of partial differential equations (PDEs) with Walsh functions offers new opportunities 
to simulate many challenging problems in mathematical physics. The approach was developed to 
better simulate hypersonic flows with shocks on unstructured grids. [1] 

The Walsh function approach is unique in that integrals and derivatives are computed using 
simple matrix multiplication of series representations of functions without the need for divided 
differences. The product of any two Walsh functions is another Walsh function - a feature that 
radically changes an algorithm for solving PDEs. 

While this algorithm holds great promise, there are coding infrastructure challenges that make 
it difficult to further explore and develop numerical simulation tools using Walsh functions. Some 
unique strengths of Walsh functions are the capabilities to systematically account for non-linear 
interactions of component waves and to explicitly propagate boundary conditions across the so- 
lution domain. However, the implementation of these capabilities presents tedious programming 
challenges. Multiplication and division of two Walsh functions are fundamentally different from any 
currently supported operations in any programming language. Keeping track of Jacobians required 
for the solution of systems of PDEs can become especially tedious in this environment. Fortunately, 
all of this tedious detail can be supported in the background through use of a module that supports 
operator overloading and other commonly used functions and routines in the solution of partial 
differential equations. This support greatly simplifies the learning curve for new users wanting to 
explore the use of Walsh functions in solving PDEs. 

The purpose of this manual is to document the Walsh function support infrastructure in FOR- 
TRAN module WALSH_TOOLS. It also provides three demonstration problems providing examples 

of how the WALSH TOOLS infrastructure is used to solve PDEs. It is assumed the reader/user has 

some familiarity with FORTRAN 2003 in general and the use of derived types in particular. This 
manual is intended to be used in association with the original documentation [1] of the algorithm. 
Only the interface to Walsh function operations and algorithms is described; the reader must refer 
to the original documentation as well as the FORTRAN program files described herein for details 
of implementation. 

Two algorithms not previously documented are supported within the module WALSH_TOOLS. 

• The Fast Walsh Transform (FWT) in multi-dimensions is key to efficient solution of problems 
using a large number ( N ) of terms in the series. An element of any solution algorithm requires 
the transformation from wave number representation to discrete value representation and back. 
The FWT implements this transformation in order Nlog(N) operations (if N = 2 P ) rather 
than order N 2 operations as documented previously. [1] 

• The Fast Walsh Reciprocal (FWR) provides the inverse of a Walsh symmetric matrix in order 
Nlog(N) operations. Any problems involving the evaluation of the reciprocal of a function 
requires evaluation of a Walsh symmetric matrix inverse. 

Documentation and derivation of the FWT and FWR algorithms are planned for an upcoming paper. 
For now, the implementation of these algorithms is documented in the accompanying FORTRAN 
code. 

The remainder of this manual is organized as follows. Section 2 documents the Walsh function 
specific oerations, functions, and subroutines that are intended to provided infrastructure support 
for new research involving Walsh functions in the solution of PDEs. Section 3 documents the 


4 


routines used to define boundary conditions and exact solutions to the demonstration problems. 
Section 4 briefly describes the linear solver software routines [2] used to solve the PDEs with Newton 
relaxation. Section 5 defines the three demonstration test problems and documents how module 
WALSH _TOOLS is used obtain their solution. Finally, Section 6 documents results of the problems 
defined in Sec. 5. The ultimate goal of this effort is to encourage new users to borrow and improve 
upon the algorithms documented here to further advance the use of Walsh functions in the solution 
of non-linear, partial differential equations. 
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2 Walsh Function Module: WALSH TOOLS 


The Walsh function module WALSH_TOOLS includes all derived types, procedures, functions, and 
subroutines required to define and use Walsh functions in the solution of partial differential equa- 
tions. It is assumed that users will want to start with the current code structure and then make 
changes associated with their own research and application needs. Only variables, functions, and 
routines declared public - intended for general use - are documented here. Note that the Jacobians 
of every variable of type walsh are automatically computed. 

2.1 Parameters and Arrays 

2.1.1 dp 

integer, parameter : : dp = selected_real_kind(15, 307) 

All real values in all modules are set to type real(dp). The parameter dp allows the programmer 
to modify the real precision throughout. Double precision is the default. 

2.1.2 gns 

type(g), dimension ( : ) , allocatable :: gns 

The array of orthonormal Walsh functions created in subroutine g_hat_bisection or subroutine 
gn_ recursion. 

2.2 Derived Types 
2.2.1 type(g) 
type g 

integer, dimension! :) , pointer 
integer, dimension!:), pointer 
integer 
end type 

Each orthonormal Walsh function g n is of type(g). This derived type includes three essential 
elements for describing the nth Walsh function. 

• The one-dimensional integer array segment is allocated n and includes the dimensionless 
segment lengths (1 or 2) corresponding to g n of Eq. 6 of Reference[l]. 

• The one-dimensional integer array factor is allocated n. The mth element of factor defines 
the mapped location k of the product g n 9m = 9k- This array corresponds to the nth column 
of the multiplication table V(n,m) of Table 2 in Reference [1]. 

• The integer level corresponds to the value of p in Eq. (3) of Referencefl]. Segment size for 
g n is determined by p in Eq. (2) of Reference[l]. 


segment 

factor 

level 
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2.2.2 type(walsh) 


type walsh 
integer 
integer 
integer 
integer 
integer 
integer 
logical 
integer, 
real (dp) , 
logical , 
logical , 
logical , 
logical , 
logical , 
real (dp) , 
real (dp) , 
real (dp) , 
real (dp) , 
real (dp) , 
real (dp) , 
end type 


dimension(4) 
dimension(4,2) 
dimension( : 
dimension( : 
dimension ( : 
dimension ( : 
dimension( : 
dimension( : 
dimension ( : , : , 
dimension ( : , : , 
dimension( : , : , 
dimension( : , : , 
dimension ( : , : , 


) , allocatable 
) , allocatable 
) , allocatable 
) , allocatable 
) , allocatable 
) , allocatable 
) , allocatable 
) , allocatable 
) , allocatable 
) , allocatable 
) , allocatable 


len 

ndep 

nbcx 

nbcy 

nbcz 

nbct 

jacb 

ns eg 

domain 

bdep 

bcx 

bey 

bez 

bet 

comp 

jac 

jacbcx 

jacbcy 

jacbcz 

jacbct 


Every dependent variable in the solution of a system of PDEs is of type walsh. This derived type 
includes all of the information required to evaluate all wave components and Jacobians across a 
domain for most unary and binary operations (+,-,*,/,**, f , and involving operands of 

type walsh. 


• The integer len equals the maximum number of segments in the domain. This value is 
inherited by any derived function of type walsh. 

• The integer ndep is the total number of primary, dependent variables in the domain. This 
value is inherited by any derived function of type walsh. If the domain is global then ndep 
is usually equal to the number of PDEs to be solved. If the domain is a boundary then ndep 
equals the total number of boundary conditions in the respective coordinate direction. 

• The integer nbcx is the total number of boundary conditions in the x (or a) direction for the 
system. This value is inherited by any derived function of type walsh. 

• The integer nbcy is the total number of boundary conditions in the y (or /3) direction for the 
system. This value is inherited by any derived function of type walsh. 

• The integer nbcz is the total number of boundary conditions in the z (or 7 ) direction for the 
system. This value is inherited by any derived function of type walsh. 

• The integer nbct is the total number of boundary conditions in the t (or t) direction for the 

system. This value is inherited by any derived function of type walsh. 
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• The logical flag jacb is true if the function is derived from (or is itself) a primary dependent 
variable or boundary condition. Often, PDEs include metric coefficients represented by Walsh 
functions across the domain that are only functions of the independent variables, in which 
case jacb is false. If a new function is derived from a Walsh function with jacb equal true, it 
too will assign jacb equal to true. 

• The one-dimensional, integer array nseg of length 4 includes the number of segments across 
the domain in the x, y, z, and t directions (or a, f3, 7 , and r directions), respectively. If a 
simulation includes only two directions then the segment numbers of the remaining directions 
are set to 1 - indicating a constant value of all functions in those directions. This array is 
inherited by any derived function of type walsh. 

• The two-dimensional, real array domain with dimension (4,2) includes the initial, constant 
value of the computational coordinate in the x, y, z, and t directions (or a, /3, 7 , and r direc- 
tions), in the (1,1), (2,1), (3,1), and (4,1) locations, respectively. It includes the terminating, 
constant value of the computational coordinate in the x, y, z, and t directions (or a, /3, 7 , 
and t directions), in the (1,2), (2,2), (3,2), and (4,2) locations, respectively. If a simulation 
includes only two directions then the initial and terminating values of the remaining directions 
are set to 0 and 1. This array is inherited by any derived function of type walsh. 

• The one-dimensional, logical array bdep has length ndep. The mth element of bdep is true 
if it is derived from the mth primary, dependent variable. Otherwise it is false. 

• The one-dimensional, logical array bcx has length nbcx. The mth element of bcx is true if 
it is derived from the mth boundary condition in the x (or a) direction. Otherwise it is false. 

• The one-dimensional, logical array bey has length nbey. The mth element of bey is true if 
it is derived from the mth boundary condition in the y ( or /3) direction. Otherwise it is false. 

• The one-dimensional, logical array bez has length nbez. The mth element of bez is true if 
it is derived from the mth boundary condition in the z ( or 7 ) direction. Otherwise it is false. 

• The one-dimensional, logical array bet has length nbet. The mth element of bet is true if it 
is derived from the mth boundary condition in the t ( or r) direction. Otherwise it is false. 

• The one-dimensional, real array comp has length len. It contains the wave components Ai 
of the Walsh function series. 

nseg(4) nseg(3) nseg(2) nseg(l) 

E E E E Aigi(x)gj(y)g k (z)g n (t ) 

n = 1 k = 1 j = 1 i = 1 

where l = i + nseg(l)(j — 1 + nseg(2)(fc — 1 + nseg(3)(n — 1))). 

• The three-dimensional, real array jac is dimensioned (len, len, ndep). Jacobian jac (l,m,n) 
contains the partial derivative of the Zth wave component with respect to the mth wave 
component of the nth primary, dependent variable. 

• The three-dimensional, real array jacbcx is dimensioned to (len, lenx, max(l,nbcx)), where: 

lenx — max(l, nseg(2) * nseg(3) * nseg(4)) 
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The element jacbcx(l,m,n) contains the partial derivative of the Ith wave component with 
respect to the mth wave component of the nth boundary condition in the x (or a ) direction. 

• The three-dimensional, real array jacbcy is dimensioned to (len,leny, max(l,nbcy)), where: 

leny = max(l, nseg(l) * nseg(3) * nseg(4)) 

The element jacbcy(l,m,n) contains the partial derivative of the lih wave component with 
respect to the mth wave component of the nth boundary condition in the y (or (5) direction. 

• The three-dimensional, real array jacbcz is dimensioned to (len,lenz, max(l,nbcz)), where: 

lenz — max(l, nseg(l) * nseg(2) * nseg(4)) 

The element jacbcz(l,m,n) contains the partial derivative of the Ith wave component with 
respect to the mth wave component of the nth boundary condition in the z (or 7 ) direction. 

• The three-dimensional, real array jacbct is dimensioned to (len,lent, max(l,nbct)), where: 

lent = max(l, nseg(l) * nseg(2) * nseg(3)) 

The element jacbct(l,m,n) contains the partial derivative of the It h wave component with 
respect to the mth wave component of the nth boundary condition in the t (or r) direction. 

2.3 Operator Overloading 

2.3.1 assignment (=) 

public : : assignment (=) 
interface assignment (=) 

module procedure walsh_assign_dd 
end interface 

A function of type walsh on the left side of the equals sign (=) is set equal to the result of the 
expression involving functions of type walsh on the right side. 

2.3.2 operator(-f-) 

public : : operator(+) 
interface operator(+) 

module procedure walsh_plus_dd 
module procedure walsh_plus_dr 
module procedure walsh_plus_rd 
end interface 

A binary operation in which a function of type walsh may be added to another function of type 
walsh or to a real number. 
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2.3.3 operator(-) 

public : : operator(-) 
interface operator(-) 

module procedure walsh_minus_d 
module procedure walsh_minus_dd 
module procedure walsh_minus_dr 
module procedure walsh_minus_rd 
end interface 

A binary operation in which a function of type walsh may be subtracted from another function of 
type walsh or from a real number. The unary operation simply takes the negative of the function. 

2.3.4 operator(*) 

public : : operator(*) 
interface operator(*) 

module procedure walsh_prod_dd 
module procedure walsh_prod_rd 
module procedure walsh_prod_dr 
end interface 

A binary operation in which a function of type walsh may be multiplied by another function of 
type walsh or by a real number. 

2.3.5 operator(**) 

public : : operator(**) 
interface operator(**) 

module procedure walsh_expo_di 
end interface 

A binary operation in which a function of type walsh may be raised to a positive integer power. 

2.3.6 operator(/) 

public : : operator(/) 
interface operator(/) 

module procedure walsh_div_rd 
module procedure walsh_div_dr 
module procedure walsh_div_dd 
end interface 

If a function of type walsh occurs on the right side of the operator /, its reciprocal is evaluated 
using a Fast Walsh Reciprocal Transform and then the reciprocal multiplies the function or real 
number on the left. If a real number occurs on the right, its reciprocal multiplies the function on 
the left. 
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2.4 Functions 

2.4.1 f from comp 


pure elemental function f _f rom_comp(x, 

type(walsh), intent(in) : 
real(dp), intent(in) : 
logical, optional, intent (in) : 
real (dp) : 


y, z, t, a, & 
filter) 
a 

x, y, z, t 

filter 

f _f rom_comp 


The value of a function a of type walsh at location (x, y, z, t) in the domain is returned as 
f from comp. If filter is present and true, the contributions from the highest family of wave 
numbers in each direction are not included in the function evaluation. This function does not utilize 
the Fast Walsh Transform (FWT) which may be more efficient in many circumstances. 


2.4.2 absw 

function absw(a) 

type(walsh), intent(in) :: a 
type (walsh) : : absw 

The absolute value of a function a of type walsh is returned. A Fast Walsh Transform (FWT) is 
used to convert from wave number component space to functional value space, the absolute value 
of the function is computed, and then another FWT transforms the result back to wave number 
component space. 


2.4.3 sqrtw 

function sqrtw (a) 

type(walsh), intent(in) :: a 
type (walsh) : : sqrtw 

The square root of a function a of type walsh is returned. A Fast Walsh Transform (FWT) is 
used to convert from wave number component space to functional value space, the square root 
of the function is computed, and then another FWT transforms the result back to wave number 
component space. 


2.4.4 gn 

pure elemental function gn(xO,xl,x,n) 


real (dp) , 

intent (in) 

: xO 

real (dp) , 

intent (in) 

: xl 

real (dp) , 

intent (in) 

: x 

integer, 

intent (in) 

: n 

real (dp) 


: gn 


The value gn of the nth orthonormal walsh function, at location x in the domain, with lower bound 
xO and upper bound xl, is returned using Eq (4) of Reference [1]. In many cases, use of an FWT 
is more efficient. 
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2.4.5 pmap 

pure elemental function pmap(nl,n2) 
integer, intent(in) :: nl , n2 
integer : : pmap 

The product of basis functions nl and n2 map to pmap. That is: g n i{x)g n 2 {x) = g pm ap(x ) for 
domain length 1. See Table 2 of Reference [1], 

2.4.6 intx 

pure function intx(f ,fa,fb,diff ) 


type (walsh) , 

intent (in) : 

: f 

type (walsh), optional, 

intent (in) : 

: fa 

type (walsh), optional, 

intent (in) : 

: fb 

logical, optional, 

intent (in) : 

: diff 

type (walsh) 


: intx 


If diff is present and true, the returned function intx of type walsh defines the partial derivative of 
the input Walsh function f with respect to x (or coordinate a). If diff is absent or false, the returned 
function intx of type walsh defines the integral of the input Walsh function f as a function of x 
(or coordinate a). Evaluation of the derivative and integral differ primarily by the transformation 
matrix used: T> for the derivative from Eq. 52 of Reference [1] or x T f° r the integral from Eq. 44 of 
Reference [1]. A Walsh function boundary condition variable, defining the variation of a boundary 
condition at x = domain( 1 , 1) for fa or at x = domain^ 1,2) for fb, provides additional degrees of 
freedom to satisfy PDE system boundary conditions. Only one boundary condition (fa or fb) may 
be used in the function call. 

2.4.7 inty 

pure function inty(f ,fa,fb,diff ) 


type (walsh) , 

intent (in) : 

: f 

type (walsh), optional, 

intent (in) : 

: fa 

type (walsh), optional, 

intent (in) : 

: fb 

logical, optional, 

intent (in) : 

: diff 

type (walsh) 


: inty 


Same as intx but in the y (or coordinate /3) direction. 

2.4.8 intz 

pure function intz(f ,fa,fb,diff ) 


type (walsh) , 

intent (in) : 

: f 

type (walsh), optional, 

intent (in) : 

: fa 

type (walsh), optional, 

intent (in) : 

: fb 

logical, optional, 

intent (in) : 

: diff 

type (walsh) 


: intz 


Same as intx but in the c (or coordinate 7 ) direction. 
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2.4.9 intt 


pure function intt(f ,fa,fb,diff ) 


type (walsh) , 

intent (in) : 

: f 

type (walsh), optional, 

intent (in) : 

: fa 

type (walsh), optional, 

intent (in) : 

: fb 

logical, optional, 

intent (in) : 

: diff 

type (walsh) 


: intt 


Same as intx but in the t (or coordinate r) direction. 

2.5 Subroutines 

2.5.1 allocate walsh 

pure subroutine allocate_walsh(a, nseg, ndep, nbcx, & 

nbcy, nbcz, nbct, jacb) 


type(walsh ), intent (inout) :: a 

integer, dimension(4) , intent(in) :: nseg 

integer, intent(in) :: ndep 

integer, intent(in) :: nbcx 

integer, intent(in) :: nbcy 

integer, intent(in) :: nbcz 

integer, intent(in) :: nbct 

logical, intent (in) :: jacb 


This routine allocates the elements of function a of derived type walsh and initializes its scalar 
elements. The allocatable vector elements of a are initialized to zero. If, on input, the allocatable 
elements of a are already allocated, the routine checks that the allocation is correct. If not correct, 
the elements are deallocated and then reallocated. The other inputs are defined previously for 
type(walsh). The domain boundaries domain are not changed by this routine. 

2.5.2 deallocate walsh 

pure subroutine deallocate_walsh(a) 
type(walsh ), intent (inout) : : a 

This routine deallocates all allocatable elements of a. 

2.5.3 initialize walsh 

subroutine initialize_walsh(a, nseg, ndep, nbcx, nbcy,& 

nbcz, nbct, domain, b, mseg, dep) 


type(walsh ), intent (out) : : a 

integer, dimension(4) , intent(in) :: nseg 

integer, dimension(4) , intent(in) :: mseg 

real(dp), dimension(4,2) , intent(in) :: domain 

real(dp), dimension( : ) , intent(in) :: b 
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integer, 

integer, 

integer, 

integer, 

integer, 

integer, optional 


intent (in) : : ndep 
intent (in) : : nbcx 
intent (in) : : nbcy 
intent (in) : : nbcz 
intent (in) : : nbct 
intent (in) : : dep 


Primary dependent variables of type walsh are initialized in this routine. 

• a: The function of type walsh to be initialized. 

• nseg, ndep, nbcx, nbcy, nbcz, nbct, domain: See definition in Sec. 2.2.2. 

• b: A one-dimensional array of real values used to compute the initial wave components of 
a using a FWT. The b array may be dimensioned to either the total number of segments 
in the entire domain or to the number of segments on any boundary of the domain. If 
b is dimensioned according to the number of segments on a domain boundary, then that 
distribution is injected into the rest of the domain. For example, in a time-dependent problem, 
the known initial condition at t = 0 (domain(4,l)) is used to initialize the function at later 
times. 

• mseg: The number of segments for each computational direction of the initializing array b. 
The array is ordered by l = i + mseg(l) (j — 1 + mseg(2)(/c — 1 + mseg(3)(n — 1))). 

• dep: An integer specifying the identity of the initialized variable used for defining the Jaco- 
bian. If a is the second of three primary dependent variables, then ndep — 3 and dep — 2. 
If a is the fourth of six boundary conditions in the x (or a) direction, then ndep — 6 and 
dep — 4. The initialization of every primary dependent variable will have identical values of 
ndep, nbcx, nbcy, nbcz, nbct that are determined by the system of PDEs being solved. 
The initialization of every boundary condition will have nbcx — nbcy — nbcz — nbct 

0, and ndep equal to the total number of boundary conditions in that direction on both the 
initial and terminal boundaries. If dep is absent or equal to 0, then the initialized function a 
is not dependent on any dependent variable or boundary condition. 


2.5.4 fast walsh gather 4d 


subroutine fast_walsh 

type (walsh ), 
type (walsh ), 
logical , 
logical , 
integer, 


_gather_4d(a, b, need_jac, & 

direction, truncate) 

intent (in ) : : a 
intent (out) : : b 
intent (in ) : : need_jac 
intent (in ) : : direction 
intent (in ) : : truncate 


This routine is the single public interface for executing Fast Walsh Transforms (FWT). Other private 
routines exist internal to the WALSH _TOOLS module for executing parts of the transformation. 

• a: Input function of type walsh to be transformed. 
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• b: Transformed output function of type walsh. 

• need jac: A logical flag indicating the need to transform the Jacobian elements of a. Trans- 
form the Jacobians if true. 

• direction: A logical flag indicating the direction of the transformation. If true, then the 
transformation goes from discrete functional values at smallest segment centers in the do- 
main to the Walsh function wave components. If false, the transformation goes from wave 
component space to discrete values. 

• truncate: An integer flag indicating the level of truncation in the transformation. 

0 - no truncation 

1 - truncate highest family of terms in the first coordinate direction 

2 - truncate highest family of terms in the second coordinate direction 

3 - truncate highest family of terms in the third coordinate direction 

4 - truncate highest family of terms in the fourth coordinate direction 

5 - truncate highest family of terms in all coordinate directions 

2.5.5 g hat bisection 

subroutine g_hat _bi sect ion (level_max) 
integer, intent(in) :: level_max 

• level max: The maximum value of p used in any of the computational coordinate directions 
where the maximum number of Walsh functions equals 2 P . 

This routine calculates gns (see Sec. 2.1.2) by the original bisection algorithm of Reference^]. 
This algorithm does not use recursion relations. It is less efficient than the recursion algorithm in 
subroutine gn_recursion) and it does not compute the multiplication table. It is retained to enable 
independent checks of the recursion algorithm. 

2.5.6 g recursion 

subrout ine g_recurs ion ( level_max) 
integer, intent(in) :: level_max 

• level max: The maximum value of p used in any of the computational coordinate directions 
where the maximum number of Walsh functions equals 2 P . 

This routine calculates gns (see Sec. 2.1.2) by the recursion algorithm of Reference[l], 

2.5.7 check walsh 

subroutine check_walsh(a) 

type(walsh), intent(in) :: a 

This routine prints information about a function a of type walsh which is often helpful in debugging 
new code. It prints all of the scalar and small array elements of the derived type and the first 8 
elements of the larger arrays in the derived type. 
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2.5.8 setup domain 


subroutine setup_domain(sa 
real(dp), intent(in) :: 

integer, intent (in) : : 

integer, intent (in) : : 


sb , p , overlap , N , ds , sap , sbp , L) 


sa, 

P 


sb 


true endpoints 
N = 2**p 


overlap 


0 = “ 1122 , 
1 = 1 - 122 , 
2 = 11-22 


integer, 
real (dp) , 
real (dp) , 
real (dp) , 


intent (out) 
intent (out) 
intent (out) 
intent (out) 


N 


ds 

sap, sbp 
L 


number of segments 
segment length 
overlap endpoints 
sbp - sap 


The solution of PDEs may employ overlap of neighboring domains or overlap of a single domain 
across a boundary. The effective boundary location is defined by Eq. 100 of Reference [1]. 

• sa, sb: The true endpoints of the domain to be solved. 

• p: The number of Walsh functions in the direction defined by (sa, sb) is 2 P . 

• overlap: An integer flag defining the extent of the overlap. 

0 - no overlap, the effective domain is exactly the true domain. 

1-1/2 segment overlap, the midpoint of the terminating, smallest segment is positioned over 
the true boundary. 

2 - full segment overlap, the terminating smallest segment is positioned across the true bound- 
ary. 

• N: Maximum number of segments in the given computational direction. 

• ds: Length of smallest segment in the effective domain. 

• sap, sbp: Endpoints of the effective domain, including overlap. 

• L: Length of the effective domain. 
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3 Boundary Module: BOUNDARY 


The module BOUNDARY includes all functions and subroutines required to define Walsh function 
boundary conditions for the demonstration problems. 

3.1 overlap 



Figure 1: Schematic of the distribution of segments in gs(x) as a function of the overlap parameter. 


The user-defined integer flag overlap is a key element of boundary condition formulation. As 
shown in Fig. 1, overlap defines the extent to which the terminating segments of the Walsh basis 
functions overlap the boundary of the physical domain. For overlap — 0, there is zero overlap. 
For overlap = 1, the midpoint of the terminating segment coincides with the boundary of the 
physical domain. The two half segments extending beyond the boundary sum to one. For overlap 
— 2, the terminating segments bound the physical domain. In this case, the two full segments 
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extending beyond the boundary sum to two. The total length of segments extending beyond the 
domain boundaries leads to the following equation for the smallest segment length in the Walsh 
series representation of the solution. 


Ax = 


X b ~ X a 


N a — overlap 

Here N a = 2 Pa is the terminating Walsh function index in the series representation, 


N a 


q(x) = qw%comp(q)g a (x) 


a=l 


(1) 

(2) 


and the symbol % indicates that the wave number component comp is an element of the Walsh 
function series type qw as defined in Sec. 2.2.2. 

Numerical tests indicate that second-order differential equations are most accurately modeled 
with overlap ^2. This overlap requires that the Walsh series representation of the dependent 
variable(s) has a value on the bounding segments equal to the given Dirichlet boundary conditions. 
A zero gradient boundary requires the values on the two terminating segments to be equal. Examples 
of the coding for these formulations are included in the demonstration test cases for Burgers Equation 
and the Riemann problem. 

The first-order linear advection equation test case uses overlap = 1 to impose a periodic 
boundary condition. In this case, the dependent variable value on the left terminating segment at 
x a is set equal to the value on the right terminating segment at 

Figure 2 shows the distribution of segments across subdomains that span the global domain. 
Subdomains are introduced to reduce the storage requirements for Jacobians. In the example of 
Fig. 2, if a single series with 16 terms spans the global domain, then a single Jacobian with 
16 2 = 256 is evaluated. (Boundary condition contributions to the Jacobian are not included in this 
simple example.) If the global domain is subdivided into 4 equal parts (shown as black, green, red, 
and blue in Fig. 2), and each of these subdomains is spanned by a series with 4 terms, then the 
discretization of the global domain is equivalent to the single 16 term representation. However, the 
Jacobian storage is a factor of 4 smaller ( 4 x 4 2 = 64) and Newton convergence is compromised 
when the global domain is partitioned in this manner. 

When subdomains are used, inter-domain boundary conditions are required. The overlap pa- 
rameter plays the same role in the inter-domain boundaries as in the global boundaries. For overlap 
— 0, neighboring subdomains (illustrated as black to green, green to red, red to blue in Fig. 2) 
have no overlap; their endpoints are coincident. For overlap —1, neighboring subdomains have one 
overlapping segment. For overlap — 2, neighboring subdomains have two overlapping segments. 
The three light black vertical lines in Fig. 2 define the inter-domain boundaries. The extent of 
overlap of the subdomains as a function of overlap is exactly the same as shown for the global 
domain boundaries in Fig. 1. The minimum segment size when tid subdomains are utilized is given 

by 


Ax 


Xb - Xg 

nD(N a — overlap) 


( 3 ) 


The advection test case uses overlap — 1 to compute inter-domain boundary conditions. The 
advection direction is from left to right. The function value computed at the far right segment 
of the black, green, and red sections is used to define the far left green, red, and blue segments, 
respectively. The Burgers equation and Riemann test cases include dissipation terms. In these 
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Figure 2: Schematic of the distribution of segments in the case where a four term Walsh series is 
implemented across each of four subdomains spanning the global domain. 


cases, information is transferred in both directions across the inter-domain boundaries. The effect 
of this data transfer is to make the functional values at the overlapping segments of neighboring 
subdomains identical. 


3.2 profile 


function prof ile (xp , demo_code , var) 


real(dp), intent(in) 
integer, intent (in) 
integer, intent (in) 
real (dp) 


xp 

demo_code 

var 

profile 
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The initial conditions for the three test problems are returned in profile as a function of xp. 

• xp: position 

• demo code: test problem identification 

0 - advection problem 

1 - Burgers equation 

2 - Riemann problem 

• var: returned primary dependent variable 

1 - density, p for demo code = 2 , velocity, u for demo code — 0 or 1 

2 - momentum, pu for demo code = 2 

3 - energy density, pE for demo code — 2 

3.3 profile burgers 

function prof ile_burgers(x,nu) 
real(dp), intent(in) :: x, nu 
real (dp) : : prof ile_burgers 

The analytic, steady solution for Burgers equation is returned in profile burgers computed as a 
function of x and nu. 

3.4 profile_sod 

subroutine prof ile_sod(x,t,p, rho , u, e) 
real(dp), intent(in ) : : x, t 
real(dp), intent(out) :: p,rho,u,e 

The analytic, time dependent solution [3] for the Riemann problem is returned in profile sod 
computed for the Sod test conditions [4], 

• x: position, —1 < x < 1 

• t: time 

• p: pressure 

• rho: density 

• u: velocity 

• e: energy 

3.5 boundary from interior 

public : : boundary_f rom_interior 
interface boundary_f rom_interior 

module procedure boundary_f rom_interior_l 
module procedure boundary_f rom_interior_2 
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end interface 

subroutine boundary _f rom_ inter i or _1 (qwp , qb , f ace , overlap , qbj ) 
type(walsh), intent(in) :: qwp 

real(dp), dimension( : ) , intent (inout) :: qb 
integer, intent(in) : : face 

! l=xmin 2=xmax 
! 3=ymin 4=ymax 
! 5=zmin 6=zmax 
! 7=tmin 8=tmax 

integer, intent (in) : : overlap 

! 0=~1122 
! 1 = 1-122 
! 2 = 11-22 

real(dp), dimension( : , : ) , optional, intent (inout) :: qbj 

subrout ine boundary _f rom_interior_2 (qwp , qb , f ace , overlap , qbj ) 
type(walsh), intent(in) :: qwp 

real(dp), dimension( : ) , intent (inout) :: qb 
integer, intent(in) : : face 

! l=xmin 2=xmax 
! 3=ymin 4=ymax 
! 5=zmin 6=zmax 
! 7=tmin 8=tmax 

integer, intent (in) : : overlap 

! 0=~1122 
! 1 = 1-122 
! 2 = 11-22 

real(dp), dimension( :,:,:) , optional, intent (inout) :: qbj 

Boundary conditions are usually defined as a function of discrete values in physical space. Conse- 
quently, a FWT is used to transform components of primary dependent variables in sequency space 
to discrete values in physical space. A Walsh series dependent variable is denoted qw. The FWT of 
qw is denoted qwp. Both functions are of type walsh. The element comp of qwp is an array of 
discrete values at segment midpoints in the domain. The element jac of qwp is an array containing 
the Jacobian of discrete values of the dependent variable with respect to the wave components of 
the primary dependent variables. 

• qwp: The FWT of qw in the domain. 

• qb: The interpolated value of the discrete function qwp%comp on the boundary. 

• face: Boundary code. 

1 = xmin, 2 = xmax 
3 = ymin, 4 = ymax 
5 = zmin, 6 = zmax 

• overlap: An integer flag defining the extent of the overlap. 

0 - no overlap, the effective domain is exactly the true domain. 
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1-1/2 segment overlap, the midpoint of the terminating, smallest segment is positioned over 
the true boundary. 

2 - full segment overlap, the terminating smallest segment is positioned across the true bound- 
ary 

• qbj: The Jacobian of qb with respect to the wave components of qw. The argument is 
optional and if it is not present then a Jacobian does not need to be returned. In general, 
the interpolated value (qb) of the discrete function on the boundary is a function of qwp 
at two terminating segments for overlap -2. However, qb is a function of qwp at a single 
terminating segment for overlap — 0 or 1. If the entry is a two-dimensional, real array, then 
the Jacobian at the single terminating segment for overlap — 0 or 1 is returned. If the entry 
is a three-dimensional, real array, then the Jacobian at the terminating segment (third index 
= 1) and its neighbor (third index = 2) for overlap — 2 is returned. 
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4 LAPACK Module: LAPACK 

The linear solver sgesv from LAPACK [2] is required. That routine, and lower level routines required 
by sgesv have been assembled into a module to simplify installation. All real values used in these 
routines have been converted to real(dp) (see Sec. 2.1.1). Because the code is changed, the routine 
is renamed to sgesv mod per request of the authors. [2] 

LAPACK is available from netlib via anonymous ftp and the World Wide Web at 
http://www.netlib.org/lapack, and is governed by the following license: 

Copyright (c) 1992-2013 The University of Tennessee and The University of Ten- 
nessee Research Foundation. All rights reserved. Copyright (c) 2000-2013 The Univer- 
sity of California Berkeley. All rights reserved. Copyright (c) 2006-2013 The University 
of Colorado Denver. All rights reserved. 

COPYRIGHT 

Additional copyrights may follow 

HEADER 

Redistribution and use in source and binary forms, with or without modification, are 
permitted provided that the following conditions are met: 

- Redistributions of source code must retain the above copyright notice, this list of 
conditions and the following disclaimer. 

- Redistributions in binary form must reproduce the above copyright notice, this list 
of conditions and the following disclaimer listed in this license in the documentation 
and/or other materials provided with the distribution. 

- Neither the name of the copyright holders nor the names of its contributors may be 
used to endorse or promote products derived from this software without specific prior 
written permission. 

The copyright holders provide no reassurances that the source code provided does not 
infringe any patent, copyright, or any other intellectual property rights of third parties. 

The copyright holders disclaim any liability to any recipient for claims brought against 
recipient by any third party for infringement of that parties intellectual property rights. 

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CON- 
TRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUD- 
ING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABIL- 
ITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO 
EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE 
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CON- 
SEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCURE- 
MENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROF- 
ITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THE- 
ORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF 
THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF 
SUCH DAMAGE. 
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5 Demonstration Program: DEMO 

The demonstration program DEMO includes three test cases providing guidance on the use of the 
Walsh function module WALSH _TOOLS to solve partial differential equations. All of the sample 
problems are time dependent. Two are non-linear. One involves a system of equations. 

If the temporal component is computed with p T = 0, then the time derivative is approximated 
with a divided difference. That is ~ (q(-M+Ar) q(i.t)) . j n this case of a single segment spanning 
time step dt, there are not enough degrees of freedom to use the differentiating matrix V. [1] When 
p T > 0, the differentiating matrix is employed and divided differences are not needed. Sample code 
blocks will be provided to illustrate this point. 


5.1 Test Cases 

5.1.1 Advection 


The linear, first-order advection equation is 


du du 

~dt +C di^ 


(4) 


where c = 1 is the wave speed and u(x, 0) = ito(x) is the initial condition on the domain 0 < x < 1. 
The initial profile moves to the right. A periodic boundary condition is applied so that as the 
profile exits the right boundary it reemerges from the left. The test is designed to provide an 
easily evaluated metric to measure how well the initial profile is preserved. The initial profile in the 
demonstration case (see Fig. 3) is a sawtooth defined by 


Or) 


< 1 — 4\x — I 

1 


for 0 < x < \ 
for \ < x < | 
for | < x < 1 


Code for Eq. 4 using WALSH _TOOLS is written: 
if(N_tau > l)then 

resw(l) = intt (qlw(m) ,f a=ql0w(m) , dif f = . true . ) 

+ wave_speed*intx(qlw(m) ,fa=qlaw(m) , dif f=. true. ) 

else 

resw(l) = (qlw(m) - ql0w(m))/dt 

+ wave_speed*intx(qlw(m) ,fa=qlaw(m) ,diff=.true. ) 

end if 


(5) 


In this example, if N T > 1 ( p T > 0), then the time derivative term is calculated inside the function 
call to intt and it uses the differentiating matrix. If N T = 1 ( p T = 0), then the time derivative 
term is calculated as a divided difference. In both cases, the space derivative is calculated inside the 
function call to intx. Note that for this linear, first-order equation, only a single Walsh function 
dependent variable qlaw(m) is required to satisfy a boundary condition in space and qlOw(m) 
is required to satisfy an initial condition in time. Also note that index m refers to a subdomain. If 
Pdomain = 0, then m is identically equal to 1. 
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Figure 3: The initial profile for the advection test case, uq(x). 


A subtle but critically important point is encountered when defining the dependent variable 
boundary condition qlOw(m). In general, one knows the initial condition. Why not simply repre- 
sent qlOw(m) explicitly as the Walsh transform of the known, initial condition? Different variations 
of this idea have been explored with no satisfactory solution. Numerical tests indicate that the com- 
puted qlOw(m) provides a close but inexact representation of the initial condition on the segment 
center. Allowing qlOw(m) to be part of the dependent variable set provides additional degrees of 
freedom that enables the computed solution on the global domain to exactly match the specified, 
initial conditions. This same observation applies to domain boundary conditions in space that are 
explicitly known. Treating boundary conditions as dependent variables provides additional relief at 
domain corners where spatial and temporal boundary conditions overlap. Close examination of the 
source code will reveal that boundary conditions at these corners are not overspecified; rather, the 
highest component Walsh function contribution to one of the boundary condition series is set to 
zero to reflect the fact that one less degree of freedom is required to close the system. 

5.1.2 Burgers Equation 

The non-linear, second-order Burgers equation is 

du 1 du 2 d 2 u 

dt + 2 dx U dx 2 

where the diffusivity v > 0 is input by the user. (See Sec. 5.3.) The initial condition is a linear 
function uq(x) = — x on the domain —1 < x < 1 with boundary conditions u a (t) = u( — l,t) = 1 
and Ub(t) = u(l,t) = — 1 . 

Given these initial and boundary conditions, the solution profile evolves according to the local 
value of u. At any given time, points with positive value of u move to the right and points with 


( 6 ) 
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negative value of u move to the left with speed u. This motion is dissipated as a function of diffusivity 
v. Large values of v evolve profiles that are nearly linear between the boundary values. Small values 
of v evolve profiles with an abrupt, shock-like transition from 1 to -1 at x = 0. Examples will be 
provided in Sec. 6.2. 

If the diffusivity is specified 0, then an artificial diffusion term must be included to maintain a 
stable computation while accommodating boundary conditions from both the left and right. The 
artificial dissipation is defined 


= -( dx p ) 2 


du 

dx 


( 7 ) 


where dx p (given by Eq. 3) is the smallest segment size used in the Walsh series representation of 
the solution. The use of a second-order, artificial dissipation here and in the Riemann problem that 
follows presents a simple solution for communicating information from opposite boundaries. Such 
communication was not an issue in the linear advection test case in which all waves travel from 
left to right. It is thought that a characterstic-based formulation of this problem could offer better 
accuracy as v — >• 0 and the transition from u = 1 to u = — 1 occurs over a length smaller than dx p . 
For now only the artificial dissipation is used to address this issue. 

Code for Eq. 6 using WALSH _TOOLS is written: 


if(N_tau > l)then 

resw(l) = intt (qlw(m) ,f a=qlOw(m) , dif f = . true . ) & 

+ intx(0 . 5_dp*qlw(m)**2-tauxw,f a=qlbw(m) , dif f = . true . ) 

else 

resw(l) = (qlw(m) - qlOw(m))/dt & 

+ intx(0 . 5_dp*qlw(m)**2-tauxw,f a=qlbw(m) , dif f = . true . ) 

end if 


Representation of the time derivative is exactly the same as discussed above for the advection 
equation. The shear term, tauxw — vdu/dx, was computed just above this block as: 

uxw = intx(qlw(m) ,fa=qlaw(m) , dif f=. true . ) 
if (nu==0 . _dp) then 

dx2 = 0 . 5_dp*dx_ef f **2 
tauxw = dx2*absw(uxw) *uxw 
else 

tauxw = nu*uxw 
end if 


Note that if v = 0, then the artificial dissipation is added as defined in Eq. 7. Also note that this 
second-order equation has two boundary conditions that must be engaged. The Walsh function 
dependent variables qlaw(m) used in the definition of uxw — du/dx and qlbw(m) used in the 
definition of <9(.5 u 2 — vdu/ dx) / dx, provide additional degrees of freedom required to satisfy the 
boundary conditions. The Walsh function dependent variable qlOw(m) used in the definition of 
du/dt provides the additional degrees of freedom to satisfy the initial condition. 


5.1.3 Riemann Problem 

The Riemann problem in one-dimension defines the compressible gas dynamic flow that follows 
the breaking of a virtual diaphragm separating two, constant initial states. The solution typically 
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involves the propagation of a shock, a contact discontinuity, and an expansion fan. A particular 
instance of the Riemann problem was defined by Sod [4] and is frequently used in the validation of 
computational fluid dynamics codes. 

The one-dimensional, time-dependent, compressible gas dynamic equations are: 


dpu 

~dt 


+ 


dp dpu 
dt ~l~ dx 
d(p + pu 2 ) 


dpE 

~dt~ 


+ 


dx 

dpuH 

dx 


P 


e 

H 


0 


( 7-1 )pe 

E - — 

2 

E + P - 


( 8 ) 

( 9 ) 

(10) 

( 11 ) 

(12) 

(13) 


The need for artificial dissipation was discussed in the previous section. It is added in Eqs 14 - 

16. 


dpu 

~dt 


dp 

dt 

d_ 

dx 


+ 


d 

dx 


pu — zq 


p + pu - V 2 


dp 

dx 

dpu 

dx 


dpE d 
dt dx 


puH — Zq 


dpE 

dx 


0 

0 

0 


The artificial diffusion coefficients defined in Eqs. 17 - 19 are of order ( dx p ) 2 . 


(14) 

(15) 

(16) 




V2 


"3 


(< dxp) 2 
(' dx p ) 2 
(■ dxp ) 2 


dp 


dx 

dpu 


dx 

dpE 


dx 


+ 0 dx p ) 2 
+ ( dx p ) 2 
+ ( dx p ) 2 


(17) 

(18) 
(19) 


The equations are solved on the domain — 1 < x < 1 with 7 = 1.4. Constant initial conditions 
on the left (x < 0) and right ( x > 0) for the Sod test case [4] are given by 


PL = 1 
PL = 1 

ul = 0 


Boundary conditions are given by 

% 

dx 

pu(— 1 , t) = 0 
dpE 


(-M) = 0 


dx 


(-M) = 0 


PR = 0.1 
PR = 0.125 
UR = 0 


dp 


(M) = 0 


dx 

pu( 1 , t) = 0 
dpE 


dx 


(M) = 0 
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Code for Eq. 14 - 16 using WALSH_TOOLS is written: 


if(N_tau > l)then 

resw(l) = intt (qlw(m) ,f a=qlOw(m) , dif f = . true . ) & 

+ intx(q2w(m) -qlxw, f a=qlaw(m) ,diff=.true. ) 
resw(2) = intt (q2w(m) ,fa=q20w(m) , dif f= . true . ) & 

+ intx(pw + rhou2w - q2xw,fa=q2aw(m) , dif f= . true . ) 
resw(3) = intt (q3w(m) ,fa=q30w(m) , dif f= . true . ) & 

+ intx(q2w(m) *htw - q3xw,fa=q3aw(m) , dif f=. true. ) 

else 

resw(l) = (qlw(m) - qlOw(m))/dt & 

+ intx(q2w(m) -qlxw, f a=qlaw(m) ,diff=.true. ) 
resw(2) = (q2w(m) - q20w(m))/dt & 

+ intx(pw + rhou2w - q2xw ,fa=q2aw(m) ,diff=.true. ) 
resw(3) = (q3w(m) - q30w(m))/dt & 

+ intx(q2w(m) *htw - q3xw,fa=q3aw(m) , dif f=. true. ) 

end if 


The main difference between this code sample and previous ones in this section is that there are 
three dependent variables and three corresponding residual equations to be solved. 

5.2 Compilation 

Any modern FORTRAN compiler should be satisfactory for creating the executable. In the present 
example, the FORTRAN compiler “gfortran” is used. Create the following four files in a working 
directory and enter the command “make” to create the executable “demo”. 

5.2.1 Makefile. env 

SHELL = /bin/sh 
F90FLAGS = -02 -g 
F90C0MPILER = gfortran 
F90LINKER = gfortran 
F90LIBS = -lm 

5.2.2 make, dependencies 

demo.o: demo. f 90 \ 
walsh_tools . o \ 
lapack.o \ 
boundary . o 

walsh_tools . o : walsh_tools . f 90 
lapack.o: lapack.f90 \ 
walsh_tools . o 

boundary. o: boundary. f 90 \ 
walsh_tools . o 
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5.2.3 Makefile. rules 


Note that the lines with eight leading spaces is a “tab” key entry. 

FLAGSFILE= . compileFlags 
## object file rules 
.SUFFIXES : 

.SUFFIXES : .c .F90 ,f90 .o 
. f 90 . o : 

$ (F90C0MPILER) -c $(F90FLAGS) $*.f90 

. F90 . o : 

Oecho "$(PKGFLAGS) " > $ (FLAGSFILE) 

$ (F90C0MPILER) -c $(F90FLAGS) $ (PKGFLAGS) $*.F90 

. c . o : 

$ (CCOMPILER) -c $ (CFLAGS) $ (PKGFLAGS) $ (PKGINCLUDE) $* . c 
5.2.4 Makefile 

Note that the lines with eight leading spaces is a “tab” key entry. 

include Makefile. env 
include Makefile. rules 
TARGET=demo 
default : 

$ (MAKE) $ (TARGET) 

$ (TARGET) : demo . o 

$ (F90LINKER) $(LDFLAGS) -o $0 *.o $(F90LIBS) 

clean: 

-rm -rf *.o *.stb *.mod $ (TARGET) 
include make . dependencies 

5.3 Inputs 

Upon entering the command “demo” the following prompts are provided. 

5.3.1 demo code 

Enter code for demo. 0=Advection, l=Burgers, 2=Riemann 

The answer to this prompt is stored as variable demo code and it controls the selection of the 
demonstration problem. 

5.3.2 nu 

If demo code = 1 the next prompt is 

Enter value for diffusivity (0 for inviscid) 

A value of v > 0 is required. If 0 < v < (dx p ) 2 , then the expected shock transition thickness is 
less than the smallest segment size and the solution may not converge. In this case, the artificial 
diffusivity z/q as defined in Eq. 7 is engaged by specifying v = 0. 
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5.3.3 p alpha 

The next prompt is 

Enter power of 2 for series g_alpha(x) 

An integer value for p a is required. The number of terms in the Walsh series solution in the x 
direction is N a = 2 Pa . 

5.3.4 p tau 

The next prompt is 

Enter power of 2 for series g_tau(t) 

An integer value for p T is required. The number of terms in the Walsh series solution in the t 
direction is N T = 2 Pt . 

5.3.5 p domain 

The next prompt is 

Enter power of 2 for number of domains spanning x 

An integer value for Pdomain is required. The number of subdomains across the x direction is 
rif) = 2 Pdomain as discussed in Sec. 3.1. If Pdomain = 0 the subdomain is identical to the global 
domain. 

5.3.6 overlap x 

The next prompt is 

Enter code for overlap of x-domains: 0=~1122, 1=1~122, 2=11~22 

An integer value for overlap x is required. Guidance for selection of integer code 0, 1, or 2 has 
been provided in Sec. 3.1. 

5.3.7 overlap t 

The next prompt is 

Enter code for overlap of t-domains: 0=~1122, 1=1~122, 2=11~22 

An integer value for overlap t is required. Guidance for selection of integer code 0, 1, or 2 has 
been provided in Sec. 3.1. 

5.3.8 dt 

The next prompt is 
Enter timestep 

A real value for the time step dt is required. This entry defines the domain size in t. 
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5.3.9 t max 


The next prompt is 
Enter total time 

A real value for the total time to be simulated t max is required. The solution is advanced in 
time with time step dt until total time t max is attained. 

5.3.10 truncate 

The next prompt is 
Enter truncate: l=yes, 0 = no 

An integer value for the truncation flag truncate is required. Truncation, if selected, truncates 
the contribution of the highest family of Walsh functions after each time step, damping solution 
oscillations in the vicinity of shocks. 


5.4 OutputFiles 
5.4.1 Screen Output 

The error norm (Eq. 20 ) is printed as a function of time. The norm is computed relative to the time 
accurate solution for the Advection equation and the Riemann problem. It is computed relative to 
the steady state solution (t » 0) for Burgers equation. Note that variable names with subscript 
"e" represent exact solutions to the test problems. 


Eq=i | u e (x a ,t T ) - u(x a , t T )\dx p 
Ea=l \u e (x a ,Oo) - u(x a ,t T )\dx p 


error norm = 


[\Pe(x on t T ') p[x a , t T )\ / p e (x a , t T ) 

+ \p e (x a ,t T ) ~ p(x a ,t T )\/p e (x a ,t T ) 
+ \Ue (x a , t T ~) u(x a , t T ') |] dXp 


for the Advection test case 
for the Burgers test case 


for the Riemann test case 


(20) 


The LI norm metric for convergence of wave component magnitudes and integration constants 
is printed as 11 norm for every Newton relaxation step. It is the sum of the absolute value of 
the change in these dependent variables for each relaxation step divided by the total number of 
dependent variables. 


5.4.2 Plot Files 

The solution is output for plotting using ASCII format Tecplot [5] hies. A constant over segment 
(COS) format is used in some hies to emphasize the Walsh function provenance of the solution. This 
format produces figures that appear as a stair step distribution where the width of the stairs equals 
dXp. Other hies use a segment centered (SC) format in which the exact solution at the segment 
center is compared to the Walsh function solution at the segment center. 
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• advection.dat: x, u in COS format for the advection test case. 

• advection_exact.dat: x, u e , u in SC format for the advection test case. 

• burgers.dat: x, u in COS format for the Burgers Equation test case. 

• burgers_exact.dat: x, u e , u in SC format for the Burgers Equation test case. 

• tube.dat: x, p, u, e, p in COS format for the Riemann problem. 

• tube_exact.dat: x, p e , p, u e , u, e e , e, p e , p in SC format for the Riemann problem. 

A separate Tecplot zone is defined for each time step. Animations of the solution as a function 
of time can be created by looping over zones. In the case of the Riemann problem, if the number 
of segments in the domain is greater than 32, then only every tenth time step solution is saved in 
order to reduce hie size. Note that the error norm returns a value of zero for time steps that are 
not saved. 

5.5 Flowchart 

The logical how through program DEMO is designed to enable easy contrast of function and subrou- 
tine calls for each of the three test problems. The logical how through the code is relatively simple. 
Comments in the body of the code correspond to the following list which highlight key steps in the 
algorithm. Differences between the three demonstration code cases in each section should provide 
better understanding regarding the way that WALSH_TOOLS is used. 

• Begin input and setup Parameters as function of overlap 

• Open hies for tecplot 

• Begin allocate and dehne initial conditions for each dependent variable 

• Determine maximum exponent of 2 and compute Walsh functions accordingly 

• Determine dimensions of working arrays, allocate and initialize dependent var 

• Start marching in time 

• Capture initial condition from previous time step 

• Initialize loop for Newton iterations 

• Dehne residuals for PDEs in the interior 

• Load residuals and Jacobians from the interior for linear solve 

• Gather data from interior to dehne dependent variables on the boundaries 

• Dehne dependent variables and their Jacobians on left and right boundaries 

• Dehne dependent variables and their Jacobians on initial boundaries 

• Compute residuals for boundaries 
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• Solve (drdq) dq — res 

• Update dependent variables 

• Converged or maximum number of allowed relaxations for this time step - Record solution for 
post- processing 

• Proceed to next time step if elapsed time < t_max 
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6 Sample Solutions 


6.1 Advection 

The first advection test case uses the following input: 

0 ! demo_code 

8 ! p_alpha 

2 ! p_tau 

0 ! p_domain 

1 ! overlap_x 

1 ! overlap_t 

.01 ! dt 

1 . ! t_max 

0 ! truncate 


The screen output associated with the first two time steps is: 


After 

After 

At time = 

After 

After 

At time = 


1 global relaxation steps, llnorm 

2 global relaxation steps, llnorm 

1 . 0000000000000000E-002 error norm = 

1 global relaxation steps, llnorm 

2 global relaxation steps, llnorm 

2 . 0000000000000000E-002 error norm = 


5 . 661 7743259 134005E-005 
4 . 353469 1685466348E-0 18 
4 . 827481 926308 1860E- 005 
1 . 0070258277040343E-004 
4 . 77642 11618465652E-0 18 
6 . 2407283720 14006 IE- 005 


The convergence criteria requires llnorm < 10 -10 before advancing to the next time. After 100 
time steps, completing a single cycle of the initial sawtooth profile, the output associated with the 
last two time steps is: 


At time = 

After 

After 

At time = 

After 

After 

At time = 


0.98000000000000065 error norm = 

1 global relaxation steps, llnorm = 

2 global relaxation steps, llnorm = 

0.99000000000000066 error norm = 

1 global relaxation steps, llnorm = 

2 global relaxation steps, llnorm = 

1.0000000000000007 error norm = 


1 . 633407 1728664030E-003 
1 .2065793876161089E-004 
4 . 4603672023 197547E-0 18 
1 . 66093757 12891763E-003 
9 . 7874071929656260E-005 
4 . 028867009582808 1E-0 18 
1 . 659565 18 159 10024E- 003 


The Walsh function solution of the advected profile after one cycle is compared to the exact 
solution in Fig. 4a. If a Courant number is defined as the ratio of the distance traveled by a wave 
in time step dt divided by the smallest segment size dx p , then the Courant number in this case is 
cdt/dxp = 2.55. The error norm after one cycle equals 1.66 10~ 3 . Some oscillation is evident at the 
front foot of the profile. If the case is rerun with truncate = 1 to smooth oscillations, the error 
norm after 1 cycle increases to 1.37 ICC 2 and the profile is shown in Fig. 4b. Truncation eliminates 
the oscillations at the expense of smoothing out the discontinuities at the tip and feet of the profile. 

If temporal degrees of freedom are exchanged for spatial degrees of freedom by decreasing p T to 
0 and increasing p a to 10, then the Courant number increases to 10.23 and the profile appears in 
Fig. 4c. Note that p T = 0 requires overlap t — 0. The temporal truncation error overwhelms any 
benefit from a factor 4 finer resolution in x in this case. 
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(c) p a = 10, p T = 0, truncate = 0, dt = 0.01 (d) p a = 6, p T = 6, truncate = 0, dt = 1. 

Figure 4: Advection test problem showing profile after one complete cycle at t = 1. Black is the 
Walsh function series solution. Red is the exact solution. 


In stark contrast to the previous case, if spatial resolution is sacrificed for increased temporal 
revolution by setting p a = 6 and p T = 6 and increasing the time step dt by a factor of 100, then a 
complete cycle is computed in a single time step as shown in Fig. 4d with zero error norm for the 
linear problem. Zero error means that the advected profile exactly crosses the segment midpoints. 
The Courant number in this case is 63, using the previous definition. A more representative defini- 
tion of the Courant number is to consider the distance traveled by a wave in the smallest temporal 
segment size c(dt) p divided by the smallest segment size in space ( dx ) p , in which case the Courant 
number is equal to 1. This result exhibits a resonance in that the computed profile exactly repeats 
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every time step, so that the llnorm — 0 in a single relaxation step, because the initial condition 
at the beginning of the time step equals the converged condition at the end of the time step. This 
convergence behavior is captured in the following screen output for this case. 


After 

After 

At time = 

After 

At time = 

After 

At time = 


1 global relaxation steps, llnorm 

2 global relaxation steps, llnorm 

1.0000000000000000 error norm = 

1 global relaxation steps, llnorm 

2.0000000000000000 error norm = 

1 global relaxation steps, llnorm 

3.0000000000000000 error norm = 


9 . 00049479071 15784E-004 
2 . 8991809556234448E-018 
6 . 3363534456 1093 10E- 017 
2 . 7318494984685376E-018 
1 . 7762526354977382E-016 
2 .498 11901 595673 15E-0 18 
3. 215201 1681236220E-016 


6.2 Burgers Equation 

The first advection test case uses the following input: 


1 

.1 

6 

0 

0 

2 

0 

0.1 

10 . 

0 


demo_code 

nu 

p_alpha 

p_tau 

p_domain 

overlap_x 

overlap_t 

dt 

t_max 

truncate 


The screen output associated with the first two time steps is: 


After 

After 

After 

At time = 

After 

After 

After 

At time = 


1 global relaxation steps, llnorm 

2 global relaxation steps, llnorm 

3 global relaxation steps, llnorm 

0.10000000000000001 error norm = 

1 global relaxation steps, llnorm 

2 global relaxation steps, llnorm 

3 global relaxation steps, llnorm 

0.20000000000000001 error norm = 


= 1 . 3058307 1223 1995 IE-002 

= 7 . 5467963857695639E-010 

= 6 . 3699131433843720E-018 

0.63941814680754439 
= 6 . 2273326570300310E-004 

= 4 . 6192751541354147E-010 

= 2 . 3647276144270429E-018 

0.56102491910795882 


The screen output associated with the last two time steps is: 


After 
At time = 
After 
At time = 


1 global relaxation steps, llnorm = 9 . 6066902157564443E-014 

9.9999999999999805 error norm = 5 . 8941163924331765E-004 

1 global relaxation steps, llnorm = 7 .4717043270173795E-014 

10.099999999999980 error norm = 5 . 8941163870657035E-004 


A steady solution is obtained for t > 6.6 based on attaining an llnorm < 10 -10 in a single 
relaxation step following an advance in time. The steady solution at t = 10 is shown in Fig. 5a. 
The segmented nature of the underlying Walsh function support is evident in the stair stepping 
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Table 1: Error norms for Burgers Equation with v = 0.1. 


p a error norm error norm ratio 


3 

1.28 10" 1 

- 

4 

1.33 10~ 2 

9.62 

5 

2.60 10~ 3 

5.12 

6 

5.89 10" 4 

4.41 

7 

1.40 10“ 4 

4.21 

8 

3.45 10" 5 

4.06 

9 

8.92 10~ 6 

3.87 


appearance of the solution. The stair stepping is less evident in Fig. 5b in which p a = 8 and 
segment size is a factor of 4 smaller. The error norm for this problem is recorded in Table 1. With 
each increment in p a the number of segments is doubled and the error norm decreases by a factor 
of approximately 4, indicating second-order accuracy relative to the smallest segment size. 




Figure 5: Effect of p a in case of v = 0.1, p T = 0, and Pdomain = 0 for Burgers Equation. Black is 
the Walsh function series solution. Red is the exact solution. 


As the shock thickness decreases below the width of ( dx) p with decreasing u, oscillations may 
be encountered in the solution. These oscillations are eliminated by truncating the highest family 
of Walsh function contributions to the solution at the conclusion of a time step. An example 
is presented in Fig. 6, where the solution without (a) and with (b) truncation is shown when 
v = 0.001 and p a = 8. The error norm without truncation is 7.27 10 -3 and with truncation is 
1.94 10~ 3 , indicating roughly a factor of 4 decrease in error. It is a subtle but important point to 
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note that the highest order Walsh functions still contribute to the lower-order non-linear solution 
retained after truncation through the self-mapping property under multiplication. 
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(a) truncate = 0 (b) truncate = 1 


Figure 6: Effect of truncate in case of u = 0.001, p a = 8, p T = 0, and Pdomain = 0 for Burgers 
Equation. Black is the Walsh function series solution. Red is the exact solution. 
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(b) overlap = 2 


Figure 7: Effect of overlap in case of multiple subdomains with v = 0.01, p a = 8, p T = 0, and 
Pdomain = 2. Black is the Walsh function series solution. Red is the exact solution. 
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The effect of overlap in the context of a simulation with multiple subdomains is illustrated in 
Fig. 7 for p a = 8, p T = 0, and Pdomain = 2. The interdomain boundary condition must preserve 
both the value of u and du/dx across the boundary. The overlap x = 1 environment preserves 
only the value of u at the single, shared terminating segments. The solution slope is left free to 
abruptly rotate across the shared boundary. The overlap x — 2 environment preserves the value 
of u at the two terminating segments of neighboring subdomains. By pinning two values across the 
boundary, the solution slopes are also forced to agree. Consequently, with overlap x — 1, the 
solution at the subdomain boundaries essentially stays frozen at the inititial condition of the left 
boundary for waves moving to the right {u > 0), and stays frozen at the initial condition of the 
right boundary for waves moving to the left {u < 0). The initial condition was a linear variation 
such that u(x, 0) = uq(x) = —x. (See Sec. 5.1.2.) 

6.3 Riemann Problem 

The Riemann test case uses the following input: 

2 ! demo_code 

7 ! p_alpha 

0 ! p_tau 

3 ! p_domain 

2 ! overlap_x 

0 ! overlap_t 

0.001 ! dt 

2. ! t_max 
0 ! truncate 


The screen output for the first time step is: 


After 

1 

global 

relaxation 

steps , 

llnorm 

= 

2 . 8007668595445945E-003 

After 

2 

global 

relaxation 

steps , 

llnorm 

= 

3 . 2658025447789479E-004 

After 

3 

global 

relaxation 

steps , 

llnorm 

= 

1 .4103432439592479E-004 

After 

4 

global 

relaxation 

steps , 

llnorm 

= 

1 .6225866477887173E-005 

After 

5 

global 

relaxation 

steps , 

llnorm 

= 

5 . 4659371 6749586 19E-006 

After 

6 

global 

relaxation 

steps , 

llnorm 

= 

1 . 02973 14400298979E-006 

After 

7 

global 

relaxation 

steps , 

llnorm 

= 

2 . 35222051497 19938E-007 

After 

8 

global 

relaxation 

steps , 

llnorm 

= 

5 . 2232478 176224048E-008 

After 

9 

global 

relaxation 

steps , 

llnorm 

= 

1 . 3932483790545 194E-008 

After 

10 

global 

relaxation 

steps , 

llnorm 

= 

1 .8865180918122111E-009 

After 

11 

global 

relaxation 

steps , 

llnorm 

= 

4 . 3475882498435346E-010 

After 

12 

global 

relaxation 

steps , 

llnorm 

= 

1 .7182955939989943E-010 

After 

13 

global 

relaxation 

steps , 

llnorm 

= 

1 .0161033064457939E-010 

After 

14 

global 

relaxation 

steps , 

llnorm 

= 

6 . 8963576126975439E-011 

The solution for density at t 

= 0.42 is presented 

in Fig. 8. 

The expansion moving to the left 


is evident for —.5 < x < 0. The contact discontinuity at x ~ 0.35 and the shock at x ~ 0.7 move 
to the right. The constant states between the shock and the contact discontinuity, and between 
the contact discontinuity and the head of the expansion, are in excellent agreement with the exact 
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solution. The tail of the expansion at x ~ —0.5 is slightly rounded. The contact discontinuity 
appears more dissipated than the shock. The largest differences are thought to be associated with 
the direction of characteristics approaching the discontinuities. 



Figure 8: Density profile at t = 0.42 for Riemann problem. Black is the Walsh function series 
solution. Red is the exact solution. 
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7 Summary 


The module WALSH _TOOLS has been developed and documented here to enable new users to explore 
and advance orthonormal Walsh function analyses of nonlinear, partial differential equations. The 
module contains a fundamental set of operators, functions, and subroutines that enable a new user 
to operate with Walsh functions in much the same way that Fourier analysis tools are currently 
utilized. The major difference of the Walsh function approach for nonlinear problems is the property 
of closure under multiplication - the product of any two Walsh functions is exactly another Walsh 
function. 

Three test problems documented herein include linear advection, Burgers Equation, and a Rie- 
mann problem using Sod test conditions. The first case provides examples of a profile with discon- 
tinuities in slope that is advected across the domain. The other two cases explore shock capturing 
and ability to accurately track a moving discontinuity. Code for these options appear together in 
the program so that new users may study the source code, along with narrative provided here, to 
better understand how the infrastructure in WALSH _TOOLS may be used in new applications. Of 
special note is that the infrastructure automatically provides Jacobians needed to obtain Newton 
relaxation of all equations and boundary conditions across the domain. 

The algorithms described herein should be considered a starting point for applying Walsh func- 
tion based simulations of nonlinear, partial differential equations. Improvements are expected and 
will be reported as they become available. 
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